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Extra-galactic high-energy transients: event rate densities and luminosity 

functions 

Hui Suni’2, Bing Zhang^’^.s^ zhuo Li^’^ 

ABSTRACT 

Several types of extra-galactic high-energy transients have been discovered, which include 
high-luminosity and low-luminosity long-duration gamma-ray bursts (GRBs), short-duration 
GRBs, supernova shock breakouts (SBOs), and tidal disruption events (TDEs) without or with 
an associated relativistic jet. In this paper, we apply a unified method to systematically study 
the redshift-dependent event rate densities and the global luminosity functions (ignoring redshift 
evolution) of these transients. We introduce some empirical formulae for the redshift-dependent 
event rate densities for different types of transients, and derive the local specific event rate density, 
which also represents its global luminosity function. Long GRBs have a large enough sample to 
reveal features in the global luminosity function, which is best characterized as a triple power law. 
All the other transients are consistent with having a single power law luminosity function. The 
total event rate density depends on the minimum luminosity, and we obtain the following values 
in units of Gpc“^ yr“^: O.SIqu for high-luminosity long GRBs above 10®° erg s“^; Ifidlg® for 
low-luminosity long GRBs above 5 x 10^° erg s“^; and S.SIq g above 10®° erg s“^ 

for short GRBs with three different merger delay models (Gaussian, log-normal, and power law); 

^ 10^ above 10'*'^ erg s“^ for SBOs, x 10^ for normal TDEs above 10'^'* erg s“^; 

and 0.03;^ !(2 above 10^® erg s ^ for TDE jets as discovered by Swift. Intriguingly, the global 
luminosity functions of different kinds of transients, which cover over 12 orders of magnitude, are 
consistent with a single power law with an index of -1.6. 

Subject headings: gamma-ray burst: general-stars: luminosity function-stars: fiare-supernovae: 
individual (SN 2006aj, SN 2008D) 


1. Introduction 


Extra-galactic high energy transients are intense cosmological transients whose electromagnetic emission 
peaks in the X-ray or y-ray bands. The study of extra-galactic high-energy transients has remained an active 
field in astrophysics over the years. Wide-field y-ray detectors dedicated to study y-ray bursts (GRBs) have 
led to discoveries of other types of high-energy transients, such as supernova shock breakouts (SB Os) and jets 


from tidal disruption events (TDEs). Upcoming wide field X-ray telescopes (e.g. Einstein Probe. lYuan et al 


2015!) are expected to significantly enlarge the sample of the known high-energy transients, and probably 


discover new types. 
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GRBs are the main extra-galactic 7 -ray transients. Their durations, usually described by T 90 , range 
from milliseconds to thousands of seconds. Thanks to the extensive observations led by a list of 7 -ray 
telescopes, such as BATSE, HETE-II, INTEGRAL, Swift and Fermi, our understanding of GRBs has been 
greatly advanced. Two main types base d on their durations are short GRBs (or SGRBs) with Tgo < 2 s, and 
long GRBs (or LGRBs) with Tgo > 2 s (|Kouveliotou et al~ 1993 1. Observations show that these two types of 


GRBs have distinct physical origin j^: massive star core collapses for LGRBs vs. compact star mergers for 
SGRBs (see Kumar fc Zhang 20151 for a recent review). Afterglow observations led to measurements of the 


redshifts of GRBs, allowing one to access the energetics of these events. Most LGRBs are found to have a 
typical isotropic luminosity (1-10^ keV) in the range of 10®^ ~ 10®^ erg s“^, which are called high-luminosity 
long GRBs (HL-LGRBs). A small fraction of the observed LGRBs, on the other hand, have been detected 
with peak luminosities less than 10^® erg s“^. Most of these events have distinct observational properties, 
such as long duration, smooth, single-pulse lightcurves, and are usually referred to as low-luminosity long 


HL-LGRBs ( 

Soderberg et al. 2006). More imoortantlv. they are found to form a distinct comoonent in the 

GRB luminosity function ( 

Liang et al.ll2007 ; 

Virgili et al.ll2009), suggesting that they have a distinct physical 

origin. Recent studies (e.g. Gamnana et al. 

2006 

: Waxman et al. 2007: Wang et al. 2007: Bromberg et al. 


2011 


Nakar fc Saril 20121) suggested that at least some LL-LGRBs may be related to breakouts of trans- 


relativistic shocks from exploding massive stars. 

Lacking a sensitive wide-held X-ray camera, the study of X-ray transients is only in its babyhood. 
Nonetheless, a few types of extra-galactic X-ray transients have been discovered. TDEs, bright X-ray/UV 
hares generated when super-massive black holes tidally disrupt stars, have been discovered from the archival 
X-ray survey data of the missions such as ROSAT, Chandra, and XMM-Newton (e.g. Bade et al. I 1994 


Komossa fc GreineilllOOOll . These transient s typically last for mo nths to years, much longer than the duration 
of GRBs. The discoyery of Sw J1644-I-57 (IBurrows et al.ll201lll by the Swift satellite suggested that some 
TDEs can haye super-Eddington luminosities, which point towards a relatiyistic jet associated with the TDE 
eyent. The discoyery of a thermal component associated with the prompt emission of the LL-LGRB 060218 
led to the suggestion that the signal may be related to an SBO. A ser endipitous discoyery o f an X-ray 
outburst (XRO) 080109 associated with a nearby supernoya SN 2008D (|Soderberg et al.l 1200811 suggested 
that there are indeed high-energy transients (eyen though with a much lower luminosity than LL-LGRBs) 
associated with SBOs. This discoyery established SBOs as a new type of extra-galactic high-energy transient. 


There are several important questions regarding these transients: What are the event rate densities of 
them (i.e. how often do they occur per unit volume)? How do these event rate densities depend on redshift? 
What are the luminosity functions of these transients? Do the luminosity functions evolve with redshift? 
Addressing these questions are essential to understand the progenitor systems of these transients and their 
cosmological evolution. A cross comparison among different transients may also shed light into possible 
common underlying physics behind these apparently different events. 


In the literature, some studies have b e en carried out to address th ese questions for individual transients 




Liange^^ 


Virgili et al 


2007; 


2011 


Virgih^^al 


Wanderman fc Pira 


2009; 

”3 


for normal TDEs; and Burrows et al. 


Wanderman fc Piran 


20 


2011 


4 for SGRBs; 


2008 

for SBOs; 

EsQuei et al. 

2008 


for jetted TDEs). Due to the small sample size of some types 
of these transients, the estimates of their event rate densities are usually subject to large uncertainties. 


^The duration classification sometimes leads to false identification of the physical category of GRBs, see 


IZhang et al.l 1I2OO9I I 


for a detailed discussion. 
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The total event rate density of a particular transient depends on the minimum luminosity and the shape 
of the luminosity function, which is usually not well constrained. Also, the detectors’ sensitivity, search 
algorithms, as well as instrumental selection effects all introduce additional uncertainties to the problem. 
The calculations of the intrinsic event rate density rely on the sensitivity, field of view, and working period 
of the detectors. Since these transients have been detected using very different detectors with different sets 
of parameters, special care needs to be taken in order to obtain robust results. 


In order to study the evolution of luminosity functions, one needs a large enough sample that cover a wide 
redshift range, with each redshift bin having enough events to construct a statistically meaningful luminosity 
function in the redshift bin. The X-ray transients we are studying mainly reside in the nearby universe, so that 
their redshift evolution, if any, cannot be investigated. We therefore mainly focus on the luminosity function 


luminosity function of long GRBs (e.g. 

Llovd-Ronninff et al. 

2002! Yonetoku et al. 

2004: Kocevski & Liane 

2006: Salvaterra et al. 2009. 2012: Petrosian et al. 

2OI5I 

Yu et al. 2OI5I: Pescalli et al. 

2015). Usine either 

a flux-limited sample or a non-parametric method 

(Efron & Petrosian 1992) to account for the truncation 


effect, these studies suggeste d that the LGRB data are consistent with having a redshift-evolving luminosity 


funct ion. Some papers (e.g. iLlovd-Ronning et al.l 120021 : IPetrosian et al.ll2015t lYu et al.ll2015l : IPescalli et al 


201511 suggested that the data are consistent with the hypothesis that the luminosity function is a broken 
power law with a universal shape (i.e. the power law indices before and after the luminosity break), but the 
break itself has a redshift evolution in the form ofLbOc(l-l-z)^, with k ^ 2.3. 

In this paper, we apply a unified method to systematically study the redshift-dependent event rate 
densities and the global luminosity functions (i.e. luminosity functions derived ignoring possible redshift 
evolution) of several known extra-galactic high-energy transients. For GRBs, thanks to their large sample 
size, we also investigate their luminosity function evolution. Compared with previous studies, our analysis 
has a larger sample for most transients (especially for LL-LGRBs, and TDEs), and more interestingly, we 
will derive the global event rate density distributions of all transients for the first time. In Section 2, we 
introduce the general methodology of calculating the event rate density and the luminosity function of any 
type of transient. We then introduce redshift distributions of various transients in Section 3, especially the 
new empirical models for short GRBs and TDEs. In section 4, we describe our data of all extra-galactic high 
energy transients. The results for individual events are presented in Section 5.1 (LGRBs), 5.2 (SGRBs) 5.3 
(SBOs), and 5.4 (TDEs), respectively. In 5.5, we present the global distribution of the luminosity functions 
of all the transients. Conclusions are drawn in Section 6 with some discussion. Throughout the paper, the 
concordance cosmological parameters presented by the Planck Collabor ation, i.e. Hq = 67.8kms“^ Mpc“^, 
flm — 0.308, Da = 0.692, are adopted dPlanck Collaboration et ^ 2015 ). 


2. Methodology 
2.1. Global Luminosity Function 

Quantifying the redshift-dependent event rate density and luminosity function of a certain type of 
transient is a challenging task. The observed events and their redshift and luminosity distributions are 
the results of the convolution of both (likely redshift-dependent) luminosity function and intrinsic redshift 
distribution of the events, with the proper correction of the instrumental sensitivity threshold, field of view, 
and operational time. All these complications may be disentangled with a large enough observational sample. 
However, for most transients discussed in this paper (except GRBs), the number of observational sample is 
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too small to perform such a task. In order to cross-comparing various types of transients, in this paper, we 
first ignore the possible redshift-evolution of the luminosity functions of all the transients and use the data to 
construct luminosity functions of each type of transients. This allows us to separate luminosity function and 
redshift distribution f{z). In principle, the luminosity function could be redshift-dependent. For example, 
for long GRBs for which we have collected a large enough sample, evidence of such an evolution effect has 
been collected. As a result, the luminosity function we construct in this paper only carries the meaning in 
the “global” sense. We hereby define all the luminosity functions constructed without considering redshift 
evolution as global luminosity functions (GLF). In order to use the data to construct the GLF, we also assume 
that the events with a same luminosity share the same other properties (e.g. spectral properties and detector 
parameters). This makes the observed events good indicators of the underlying general population. Also, 
since there is a wide range of the spectral peak energy (Apeak) distribution for different transients studied 
in this paper, we try as much as possible to apply the fc-corrected bolometric luminosity (I — 10^ keV) using 
the measured spectral properties of the transients ('Eg. (1^ 1. The only exceptions are the TDEs detected 
by ROSAT, Chandra, and XMM-Newton whose narrow bandpass does not allow a precise inference of the 
global spectral parameters. For these events, we use a uniform ^-correction parameter 1.4 lEg. (1301) 1. 


For a certain type of transient, we define the local specific event rate density (local event rate density 
per unit luminosity) as 

_ dpo , . 

Po,L - - ( ) 

where L denotes bolometric luminosity at the peak time hereafter. For a detector with flux sensitivity Fth, 
field of view fl, and operational time T, the detected number of events in the luminosity interval from L to 
(L -I- dL) is 


dN = —po^LdL / = “-i^Po,L9{L)dL, 

47 r Jq 1 + z dz 47 r 


where 


9{L) = 


f{z) dVjz) 
1 + z dz 


dz. 


( 2 ) 

(3) 


The dimensionless function f{z) describes the redshift-dependent event rate density, i.e. 


Pl{z) = po,Lf{z). 


(4) 


The redshift-dependent specific comoving volume reads (for the standard AGDM cosmology) 

dV{z) c 47 rll|^ 

dz Hq {I X zY[CLm{ 1 + zY+ CIkY/‘^ 


For a particular L, the maximum redshift 2 max(T), which defines the maximum volume inside which an event 
with luminosity L can be detected, can be defined by the sensitivity threshold Fth via 


Ath = 


L 

4:7tD^ (Zmax)^ 


( 6 ) 


where fc is a correction factor, which corrects the observed flux in the detector’s energy band ( 61 , 62 ) to a 
wide band in the rest frame (e.g. 1 — 10^ keV for GRBs, see Eg. (1^ 1. 


Technically, it is easier to evaluate numbers in the logarithmic luminosity bins. Equation ([2]) can be 
also written as 

dN = -^dlogL = ^(lnlO)po.L< 7 (E)M(logA) (7) 

a log L Att 








- 5 - 


Suppose AN events are detected in a finite logarithmic luminosity bin from logL to logL + A(logL), one 
then has 

^ 47r 1 11 AN 

Po,L g{L)LA{\ogL)' 

The luminosity function of a certain type of transient can be defined as 


N{L)dL oc $(L)dL, 


(9) 


with the integration of $(L) normalized to unity, i.e. 



$(L)dL = 1, 


( 10 ) 


where Lmin and L^ax are the minimum and maximum values of the luminosity distribution. One can define 
the local event rate density above a certain luminosity L, i.e. 


47r 1 


1 


AN AL 


lo g Lm ax 

OTlnlO<?(L) A(logL)T' 


( 11 ) 


logL 


The total local event rate density is therefore 

pL 


login 


1 


AN AL 


PO — P0.>Lmi„ — 


— 4^ 4 

^ ^ OTlnlOg(L)A(logL)T' 

log L/rnin 


( 12 ) 


which depends on Observationally Lmin is not well constrained, and one usually adopts the observed 

minimum luminosity, which is the upper limit of the true Tmin- As a result, the derived po is in principle 
only the lower limit of the true value. To be specific, throughout the paper, we always specify a minimum 
luminosity whenever an event rate density is quoted. 

With the definition of po fEa. dT^ l. the specific event rate density can be also written as 


Po,L = Po^{L). (13) 

Within the framework that the luminosity function does not depend on redshift (the approach adopted in 
this paper), the redshift-dependent event rate density can be written as 

P{z) = Pof{z), (14) 

where f{z) is the redshift evolution function, the form of which depends on the properties of the transients. 

The luminosity function of a certain type of transient can be derived by displaying the specific event 
rate density po,L c>c 4>(L) as a function of L. By separating the data into different luminosity bins, we use the 
observed numbers to map the relevant po.L, and then fit the data points by several empirical model forms. 
The simplest model is a single power law (SPL) form, i.e. 

$(L)ocL-“. (15) 

If the model does not define the data well, we introduce a smoothly connected broken power law (BPL) form 



LJOii 


$(T) oc 


ujcx.2-\ —1/ci; 


(16) 
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where ai and 0.2 are the power law indices before and after the break luminosity Lf,, and w defines the 
sharpness of the break. In more complicated cases (e.g. the GLF of long GRBs), one needs another power 
law segment to fit the data, and we introduce a triple power law (TPL) form, i.e. 


$(L) oc 




where ai, 02 , and 03 are the power law indices for three segments, L{, i and Lf ,,2 are the two break lumi¬ 
nosities, and uii and UJ 2 are the sharpness parameters at the two breaks. 


2.2. Luminosity function evolution of GRBs 


GRBs have a large enough sample to study the redshift evolution effect. We investigate the evolution 
effect for LGRBs and SGRBs separately. For each class, we first separate the observed GRBs into several 
redshift bins, and then apply our method to map the corresponding local specific event rate density using the 
GRBs in that specific redshift bin onl}@. Practically, for a redshift bin (^ 1 ,^ 2 ) around a certain redshift z, 
we change the integration limits in Eq.(2) to zi and min(z 2 , 2 :max(L)), respectively, and repeat the procedure 
laid out in Section 2.1. The derived local specific event rate density and event rate density are expressed 
as Pq l and pg >l„j respectively, denoting that they are derived in the redshift bin around z. Notice that 
Pq l still denotes the local value. By applying a proper correction with te redshift-evolution function /(z) 
(see details in Section 3), one can obtain the luminosity function in the redshift bin (zi < z < Z2). If GRB 
luminosity function does not evolve with redshift, then the results derived from different z bins should remain 
the same. To minimize the truncation effect by the flux limits of the detectors, we also use a sub-sample 
with a higher threshold (with peak photon flux larger than 1.8 ph s“^ cm“^) in the derivations, and compare 
the results with two thresholds. 


Alternatively, we also repeat the analysis under the sam e assumption adopted in some recent papers 
(e.g. Petrosian et al. 2015 : Yu et al. 2015 : Pescalli et al. 20151) . i.e. fix the shape of the luminosity function 
and assume that Lf, evolves as a power law with k ^ 2.3. We can then map the luminosity function from the 
data by correcting the luminosity L at z to the “local” value, i.e. Lq = L/(1 -|- z)^'^. Applying the method 
in Section 2.1 using Lq instead of L would lead to the “local” luminosity function. 


2.3. Correction factors 


To perform our analysis, the redshift of an event is needed. For GRBs, not all events have redshift 
measurements. In order to properly account for their event rate density, one needs to correct the derived 
values based on the z-known sample by the ratio between the total number and the z-known event number. 
For HL-LGRBs and SGRBs, this correction factor is approximately 3, and we adopt it in our derivations. 


The detected rates also depend on the detect or’s spectral window . For GRBs, BATSE observations 
suggested that the short-to-long ratio is about 1:3 (jPaciesas et al.l 119991 ). whereas Swift BAT, which has a 


■^There are some overlaps in adjacent redshift bins, since we want to include more GRBs in each bin to reach a better 
constraint of the corresponding luminosity function. 
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softer bandpass, only gets a 1:10 ratio (ISakamoto et al.ll201ll : lOin et alJ 2013[| . Since Swift can essentially 
detect all BATSE LGRBs, one needs to correct for the Swift-detected SGRBs by another factor of ^ 3 for 
the inferred SGRB event rate density. This factor has been taken into account in our derivations. 


For the other types of transients, the samples are not large enough to access whether we have missed 
some events due to the imperfect instrumental spectral window or the lack of redshift measurements. We 
therefore derive the event rate density using the observed events only. 


2.4. Instrumental parameters 


The three important instrumental parameters that are relevant to our derivations are the operation time 
(T), the field of view (fl), and the flux sensitivity threshold (Tth)- Table [T] lists the three parameters of the 
high-energy detectors used to study various high-energy transients discussed in this paper. While the first 
two parameters are straightforwardly defined, the definitions of sensitivity thresholds are non-trivial. This 
is particularly true for wide field triggering detec tors such as Swift BAT, with which most of the GRBs and 
jetted TDEs were detected (e.g. iLien et al.ll2014l for a detailed description of Swift BAT trigger algorithm). 
In this paper, we adopt an approximate threshold for each detector. For Swift detected events, we also adopt 
slightly different values for different types of events. 


Most GRBs were detected through the rate trigger algorithm by BAT, whereas some low-flux events, 
such as the LL-LGRB 060218 were detected through an image trigger algorithm. The image trigger is an 
additional trigger algorithm to accumulate photons from a source in a comparably longer time to look for 
transients that are not bright enough to make a rate trigger. 


For the rate-trigger GRB events detected by Swift/BAT, we adopt slightly different values for different 
sub-types. In order to have the redshift of a GRB measured, the burst usually needs to have a bright enough 
optical afterglow. On average, the z-known GRBs are brighter and therefore have a higher flux threshold 
than the standard rate trigger flux threshold. Based on the lowest value of the flux distributions of our 
sub-sample, we adopt Fth = 3 x 10“® ergcm“^s“^ for rate-triggered LGRBs (including HL-LGRBs and 
LL-LGRB 080517), and 10“^ ergcm“^s“^ for SGRBs, respectively. 


Taking HL-LGRBs as an example, we derive the flux threshold based on the lower end of the observed 
photon flux distribution. We adopt the photon flux 0.3 phcm“^s“^ as the thresholcH. The transformation 
from photon flux to specif ic flux (Eg. (1261)') requires the information of the spectrum. We consider a typical 


Band function spectrum ([Band et al. 


with a = —1. [3 = —2.3, Apeak = 200 keV at redshift z = 1. 


This threshold photon flux is translated to Fth = 3 x 10 


erg cm ^ s ^. 


The LL-LGRBs 060218 & 100316D and the two Swift-detected TDE jets were detected through image 
trigger. In the case of image trigger algorithm, the threshold flux depends on the trigger duration Ttd, with 
the dependence Fth oc The trigger duration of GRB 060218 was about 80 tQ, with a mean flux of 

^ 2.8 X 10“® erg cm“^ s“^ (T. Sakamoto, 2015, private communication). This roughly corresponds to a 


^Out of 250 HL-LGRBs in our sample, only two have peak photon flux below 0.3 phcrn”^ s“^. One of them (GRB 070419A) 
was detected through image trigger, and the other (GRB 060123) did not trigger BAT but was detected from the BAT survey 
data. 

'^http: //gen.gsfc.nasa.gov/other/191157.swift 






















threshold flux 


i^th = 2.5 X erg cm-^ s-\ (18) 



derive po,L of jetted TDEs. In fact, Sw J1644+57 was image-triggered four times. The 64s trigger duration 
was the relevant one at the peak luminosity. 


XRO 080109 was serendipitously discovered by Swift X-Ray Telescope (XRT). We use a count rate 
of 0.03 counts/s for XRT with which source variability can be detected for a ^ 1000 s observation. This 
corresponds to a flux threshold Fth ~ 10“^^ ergcm“^ s“^. 


Normal X-ray TDEs were usually disco vered from archiva l survey data of various X-ray missions, such as 
the PSPC All-Sky Survey (R ASS) ( Voees et al. 1999 ). the XMM-Newton Slew Survey Source Cata¬ 

logue (XMMSLl) (ISaxton et al.ll2008l) . and the Chandra ACIS archival data. The exposure-time-dependent 
flux sensitivity thresholds of these three detectors are listed in Table [1] For each TDE event, we consider 
the real exposure time to determine its Eih. 


3. Redshift distribution 

The redshift distribution parameter f{z) for each type of transients is essential to infer the local (specific) 
event rate density (§2), and different types of transients may have different f(z) functions. In this section, 
we discuss this function for different types of transients in detail. 


3.1. Long GRBs and supernova shock breakouts 


Long GRBs (both HL and LL) and SBOs are associated with the deaths of massive stars. To first 
order, their redshift distribution should track the history of star formation. There i s evidence that at high-z , 
the GRB rate may exceed what the star formation history ( SFH) predicts (e.g. Li 2008t Kistler et al. 20081: 
Qin et al. 2010 : Virgili et al. 2011bl: Robertson fc Ellij 20121) . However, for a wide redshift span, the SFH is 
a good proxy of the redshift distrib ution of LGRBs. In this paper, we adopt the rough analytical model of 
SFH derived by Yiiksel et al. ( 20081) using the observational data: 


/lgrb/sbo(^) = 


(1 + ^) 


3 . 4?7 


-f 


/i + zy'’-'" , fi + z 
5000 ) 


- 3 . 57 ]' 


(19) 


where rj = — 10. At 2 : < 4, this functio n is directly derived from the SFH inferred from the UV and far-IR 
galaxy data ( Hopkins &: Beacom 20061) . which is independent on the GRB observations. At z = 5 — 7, the 
SFH is enhanced from the galaxy-constrained SFH to compensate the observed GRB excess at high-z, which 
can be explained by the deficiency of the observed low-luminosity star forming galaxies missed in surveys 
but traced by GRBs. We adopt Ea. dTOl) to study both LGRBs and SBOs. The latter are only observed in 
the nearby universe, so that the modification at high-z does not enter the problem. 


®http: //gen.gsfc.nasa.gov/other/416135. swift 
®http: //gen.gsfc.nasa.gov/other/450258.swift 
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3.2. Short GRBs 


Unlike long GRBs, most short GRBs do not directly trace star formation history. Observations suggest 
that most short GRBs are consistent with having an origin not related to massive star deaths. The leading 


2()05 

Fong et al. 

20 ld see 

Berger 


(|2014l i for a recent review). 


In order to have a merger to occur, a compact star binary system needs to go through a long inspiral 
phase defined by energy loss of the system through gravitational wave radiation. The redshift distribution of 
SGRBs therefore needs to account for an additional time delay due to inspiral with respect to the creation 
of the compact binary system, which itself traces the SFH. The distribution of the merger delay time scale 
is unfortunately not known. Practically, one assumes some empirical forms of the merger delay time scale 
distribution models, and apply the data to derive best paramete rs for the delay model. Three types of 


merg er delay time models have been discussed in the literature (e.g. IVirgili et al.ll201ll : IWanderman fc Piran 
20141 and references therein): power law (PL) de cay model, Gaussia n delay model, and log-normal delay 
mode l. Current data support either a Gaussian (IVirgili et al.l 120111) or log-normal (IWanderman fc Piran 
2014) delay model, with the (PL) model disfavored (even though not completely ruled out) ( Virgili et al.l2011 ; 


Wanderman fc Piranll2014ll . Table [2] lists the three models with the best parameters currently constrained 
by the SGRB data. 

With the consideration of the merger delay time distribution, it is difficult to construct an analytical 
model for the redshift distribution of SGRBs. Instead, we perform a series of Monte Carlo simulations based 
on the SFH and merger delay distribution models to construct several redshift distributions that correspond 
to the three delay models with the best-fit parameters. First, we randomly generate 10,000 comp a ct sta r 
binary systems with a redshift distribution tracking the SFH following the model of lYiiksel et al.l (|2008r i. 
Next, we randomly generate the merger delay time scales of all these systems based on the three merger delay 
time scale models listed in Table [5J For each model, we derive the look back time of SGRBs by subtracting 
the merger delay time from the formation time, and transfer the lookback time to redshift. We repeat the 
process 10,000 times (each with 10,000 events simulated). By averaging the results, we are able to derive 
the average redshift distribution of the simulated samples. We fit the derived redshift distribution (for all 
three merger delay models) using multiple-power-law functions and derive an empirical expression of f{z) 
for each model. The simulated results with best fit empirical models are shown in FigU] The distributions 
are normalized to unity at the local universe (z = 0). The empirical formulae of f{z) for the three merger 
delay models are: 


For the Gaussian delay model (jVirgili et al.ll2011ll . one has 


fsGRBi^) — 


(1 + z)^-°^ + 


1 + z 
0.17 


O.STtj 


1-bz 

4.12 


- 8 . 0 ?) 


4.05 


- 20 . 5 ? 7 ' 


( 20 ) 


with 7 ] = —2, which is roughly a broken power law with redshift breaks at Zi = 0.45, Z 2 = 2.0, Z 3 = 3.0. 


For the log-normal delay model (IWanderman fc Piranll2014ll . one has 


/sgrb(^) — 


(1 + ^) 


5 . 7 ?) 


1 + Z 


1 . 3 ?) 


1 


—9.5r} 


1 + z 


- 24 . 5 ?)' 


0.36 J V 3.3 y V 3.3 

with rj = —2, which is roughly a broken power law with redshift breaks at Zi = 0.35, Z 2 = 1.5, Z 3 = 2.3. 


( 21 ) 
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For the power law model (|Wanderman fc Piraijl2014r) . one has 
/«jrb(^) = (1 + 2 )^'®’’ + 


\ —1.2?7 /-I I \ — 4.4?7 /i I \ — Hr? 

l + 2 :\ f 1 z\ f 1 


2.5 


3.8 


7.7 


( 22 ) 


with r] = —2.6. This model has a wider redshift distribution compared to the first two model (due to the 
wide range of the merger delay time). It is roughly a broken po wer law with redsh ift breaks at Zi = 0.42, 
Z 2 = 3.4, Z 3 = 11.3. The SGRB data do not favor this model (|Virgili et al.N201l[l . even though it is not 
completely ruled out. 


3.3. Tidal disruption events 


The event rate density of TDEs depends both on the number density of supermassive black holes 
(SMBHs) and the event rate of TDEs per galaxy. Since the TDE rate of a particular SMBH only depends on 
the properties of the galaxy itself (e.g. stellar density near the SMBH and the mass of SMBH), on average, 
it may be reasonable to assume that there is no redshift evolution of the event rate per galaxy. As a result, 
the /tde(-z) parameter o f TDEs is mostly determined b y the evolution o f the number density of SMBHs as a 
function of redshift (e.g. Donnarumma fc Rossi 2015IB . Shankar et al. ( 2013ll constructed the mass density 
distribution models for SMBHs and AGNs by considering their growth rate and radiation efficiency. We 
apply their model to derive /tde(-z) for TDEs. TDEs can happen only when the tidal disruption radius is 
larger than the event horizon of SMBHs, which gives an upper limit of the mass of SMBHs for TDEs: 


Mbh < 1-6 X 10*^ 


M, 

Mq 


Rq) ' 


(23) 


where M* and i?* are the mass and radius of the star that is disrupted by the SMBH, both normalized 
to the solar values. We therefore exclude SMBHs with mass exceeding IO^Mq. In the left panel of FiglS] 
we present the numerical hts to the number density redshift evolution of SMBHs with two mass ranges 
(10® — IO^Mq and 10^ — IO^Mq). Assuming a contant TDE rate per galaxy, in the right panel of Eigl^l 
we prese nt the normalized TD E redshift distribution /tde(-z) derived from the numetical data based on the 
model of IShankar et al.l (|2013f) . The best-fit emperical model reads 


/tde(^) = 


(1 + ^) 


0 . 2?7 


1-kz 

1.43 


-3.2r) 


1 + Z 
2.66 


-7.0?)' 


with 77 = —2. One can see that /tde( 2 ) continues to decrease with redshift, reaching 2/5 at z 
- 10 "® at z - 6 . 


(24) 


1 , and 


4. Data 

4.1. Gamma-ray bursts 

Our HL-GRB sample is only limited to Swift GRBs. This is because it is a uniform sample whose size 
is large enough to derive a well-constrained GLF. We collect all the z-known Swift GRBs before May 6 , 
2014 (250 HL-LGRBs and 20 SGRBs). This sample consists of more GRBs than previous work by intro¬ 
ducing a lower flux threshold, which allow us to better study the features near the low luminosity end. The 
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data are downloaded from th e Swift archival table available at http://swift.gsfc.nasa.gov/archive/grb_table/ 
( Sakamoto et al. 20081 2011 ). For all the bursts, the 1-s peak photon flux and photon index are recorded. 
For HL-LGRBs, the 1-s peak photon flux is directly derived from their 1-s peak luminosity. For SGRBs, 
since their durations are typically shorter than one second, we apply photon count rate with a 64-ms res¬ 
olution to derive the 64-ms peak luminosity. We calculate their 64-ms peak photon flux based on the 
ratio of the two peak count rates with different temporal resolutions (Gp^etos and Gp^s), be. Pe 4 = 
-Pi(Gp, 64 ms/Gp 4 s)- The 64-ms light curves are from the Swift Burst Ground-Analysis Information page 
(http://gen.gsfc.nasa.gov/swift_gnd_ana.html), and the 1-s light curves are obtained through re-binning. To 
ensure the correct match at the peak, the regrouping is such that the time interval at the peak time (64-ms 
resolution) matches the one for the 1-s peak photon flux provided at the GCN Circular archive. Only a 
handful of LL-LGRBs were detected so far. Table [3] collects the information of six LL-LGRBs studied in 
this paper, which were triggered not only by Swift (GRBs 060218, 080517 and 100316D), but also by other 
instruments as well: GRB 980425 by GRGG/BATSE, XRF 020903 by HETE-I^, and GRB 0 31203 by INTE¬ 


GRAL . The peak photon fluxes of pre-Swift LL-LGRBs are adopted from GCN archives or ISakamoto et al 
(2004). The data of Swift LL-LGRBs are also taken from the Swift Table. 


The time integrated spectral information is taken from the literature (references provided in Table 
ED, described by either a single power law (PL) with photon index r(N( E) oc E ^ ), a Band function 


characterized by peak energy Apeak and two photon spectral indices a and /3 (IB and et al.lll993li . or a power 
law function with an exponential cutoff (GPL) fit, i.e. N{E) oc A“ exp(—A/Ac). For the latter two models, 
an Apeak can be derived from the peak in the spectrum. For single PL fits to most BAT spectra (due 
to the narrowness of the BAT band), it is believed that the intrinsic spectrum still has a peak energy. With 
BAT GRBs jointly detected by other wide-band detectors such as Konus/Wind and Fermi/GBM, it was 
found that there exists a rough correlati on between the BAT-band photon index F and Anpak; if A neav is 
not much beyond th e BAT energy ban d ( Zhang et al. 2007 : Sakamoto et al. 20091: Virgili et al. 20121) . The 
latest best fit reads ( Virgili et al. 20121) 

log(Apeak) = (4.34 ± 0.475) - (1.32 ± 0.129)r^^^ (25) 

with a large scatter, where is photon index (positive value) defined in the BAT band. We apply this 

scaling to estimate Apeak for those GRBs whos e A^eak is no t direc tly measured. For a consistency check, 
we have also adopted Apeak values derived by Butler et al.l (2007) for a sub-sample of GRBs (for which 
Aneak is avai l able fr om that method). By repeating the calculations, we found that the derived LF using the 
Butler et al.l (j2007r i method is similar to the LF derived using our method. 


For a GRB with peak photon flux Pp, the peak flux can be calculated through 

_ PpQEN{E)dE 
" CN{E)dE ’ 


(26) 


wher e N{E) is the photon spectrum of a GRB, which is in the form of the standard Band function ((Band et al 
1993li 

(lookev) A<(a —/3 )Ao, 


N{E) = A 


{a-P)Eo 

lOOkeV 


a —(3 


exp(/3-a) E>{a-P)Eo. 


(27) 


■^The acronym “XRF” stands for “X-ray flashes”. They are softer version of GRBs. Observations show that XRFs and 
GRB s seem to form a c ontinuum in the observational and theoretical parameter spaces llSakamoto et al1l2004l: IZhang et al.l 
I 2 OO 4 I : [Sersier et al.ll200^ . In fact, GRB 060218 can be also called an XRF. We adopt the names of these events based on the 
convention adopted in their discovery papers. 
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Here the integration limits ( 61 , 62 ) enclose the detector’s energ y window (e.g. 15-150 keV for Swift BAT). 
For short GRBs, we use the i^peak data of 13 SGRBs derived bv lLti et al.l (120151) . For other 7 SGRBs whose 
Apeak was not directly measured, we estimate Apeak using Eg. (051) . For those GRBs whose Band function 
parameters are not directly measured, we adopt typical values as a = —1 and P = —2.3 for LGRBs and 
a = —0.5 and /3 = —2.3 for SGRBs. 

In order to derive the bolometric luminosity (1 — 10^ keV in the cosmological rest frame) from the 
observed peak flux, we perform a fc-correction 


-^p,bol 


= AirOfFp ■ k 


(28) 


where is the luminosity distance. The fc-correction parameter can be expressed as 


k = 


rl0‘^/il + z) 

ll/(l + z) 


EN{E)dE 


EN{E)dE 


(29) 


4.2. Shock breakouts 


Massive stars end th eir lives in catastrophic core collapses when they run out of fuel in the center 
(|Wooslev &: Weaveiill986l) . As a massive star undergos core collapse, an outgoing shock surges through the 
star. When the optical depth of photons trapped in the shock becomes unity, a SBO occurs, which provides 
the first electromagnetic emission from a supernova event. Before the shock breaks out the star, only 
neutrinos and gravitational waves can escape. At the moment of breakout, a short, bright flash is expected, 
which peaks in ultravoilet or X-rays depending on how compact the star is (lGolgatelll975l : iKlein fc Ghevalier 
19781 : iNakar fc Sarill2010l ). The SBO signal therefore carries direct information about the very early stage 
of core collapse and provides direct constraints on the type of progenitor. Since there is no electromagnetic 
precursor to alert such an event, detecting an SBO is challenging. In the X-ray and soft 7 -ray regime for 
which our paper focuses on, there are only two confirmed SBOs detected so far. One is GRB 060218/SN 
2006aj association system, which shows an X-ray thermal component wi th a temperature of ^ 0.17keV in 
a very long duration {Tqq = 2100s), soft GRB with a smooth lightcurve ( Gamoana et al. 20061) . The other 
is X-Ray Outburst (XRO) 08 0109/SN 2008D assoc iation system, which was serendipitously detected by 
Swift/XRT on 2008 January 9 (|Soderberg et al.ll2008l) . Since no 7 -ray counterpart was detected even though 
this outburst was in the field of view of BAT before and during the burst, a GRB connection was ruled 
out. These two observations have offered a great opportunity to study the detailed properties about the 
progenitors. The fact that several other LL-LGRBs seem to share similar p roperties to GRB 060218 makes 
some authors suggest that all LL-LGRBs may be associated with SBOs ('e.g. lWang et al.ll2007HNakar fc Sari 

2OI2I) . 


We use the two confirmed SBO events (with very different luminosities) to estimate their event rate 
densities. The data of GRB 060218/SN 2 006ai are already includ ed in the Table of LL-LGRBs. The data of 
XRO 080109/SN 2008D are collected from Soderberg et al. ( 2008 ). Its X-ray spectrum was fit by a power-law 
model with a photon index of 2.61 in the band of 0.3 — 10 keV, which was used for fc-correction. 
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4.3. Tidal disruption events 


Stellar tidal disruption by a supermassive black hole has been theoretically predicted by ReesI ( 1988 1. 
When a star approaches a supermassive black hole, a tidal disruption event (TDE) would occur if the tidal 
force becomes larger than the self-gravity of the star and if the radius when this happens is outside the black 
hole event horizon. Part of the disrupted debris falls into the black hole from an accretion disk, giving rise to 
a bright flare in UV or X-ray band, which lasts for several months to one year. The first TDE was dis covered 
in a quiescent galaxy NGC 5905 during the RASS survey ( Bade et al. 1996 : Komossa fc Badel 1999ll . which 
showed the characteristic luminosity decay law L oc expected for fail-back acc retion. Later, several more 
TDEs have been detected from RASS, namely RX J1242, RX J1624, RX J1420 ( Komossa &: Greiner 19991 : 
Grupe et al. 19991: Greiner et al. 2000 ). In the recent decade, more TDE candidates have been discovered 


by XMM-Newton and Chandra, mostly by comparising survey catalogs with the arch ival data (lEsauei et al 


2007 ; Esauei et aLll2008 Maksvm et al.l2o{5 Saxton et al.ll2012 ; Maksvm et al.llioij ), as well as by serendip¬ 


itous detections (|Lin et al.ll201ll) . Another ROSAT source was identified as a TDE candidate as it disap¬ 


peared in the subsequent observations with XMM-Newton and Chandra ( Gappelluti et al. 20091) . Right now 
about a dozen X-ray TDE candidates have been discovered. All these TDEs hav e large amplitudes and soft 
X-ray spectra, whose host galaxies show no sign of AGN activity ( Komoss^l2012l ). Their observed : 
luminosities range from 10"'’^ — lO^^ergs”^. 


maximum 


Recently two special TDEs were detected by Sw ift. These two TD Es, i.e. Swift J1644-I-57 (IBloom et al 


2011 : Burrows et al. 2011 ) and Swift J2058.4-1-0516 ( Genko et al. 20121) . showed some distinct features. The 
peak luminosity of Sw J1644-I-57 was around 10'^® ergs“^, and the event was followed by a radio counterpart 
(IZauderer et al.ll201ll) . Due to the super-Eddington natur e of the events, these TDEs have been in terpreted 
as relativistic jets launched from the central black hole (|Bloom et al.ll201ll : iBurrows et al.l 120111) . It has 
been claimed that ther e is a low probability that normal TDEs also host a jet similar to Swift J1644-I-57 
( van Velzen et al. 201.3h . However, it remains unknown whether two types of TDEs are indeed intrinsically 
different from others, and if so, what could be the main reason to make the difference. By measuring the 
spin parameter of the central black holes of the two Swift TDEs within the theoretical framework of the 
Blandford-Znajek mechanism. iLei fc Zhand (|20Ill) found that both black holes carry a moderately high spin. 
They then suggested that black hole spin may be the key factor to make the dichotomy of TDEs, and only 
black holes with rapid spin can launch relativisti c jets during TDEs. Further modeling su ggests that the 
j et model can successfully account for the X-ray (ILei et al.ll20I3t iTchekhovskov et al.N2014l ) and the radio 
( Metzger et al. 20I2I : Wang et al. 2014; Liu et al. 2015 ) data. 


Table [5] lists the data for all the TDE candidates. TDE flares last much longer than GRBs and SBOs. 
They usually have a relatively fast rising phase, reach and stay at the peak for some time, and then decay 
rough ly with a power law L oc The luminous state usually lasts for several months to one year (Rees 

19881) . In our sample, only two events may have been detected both in the rising a nd declining phases, so 
that the peak luminosity was measured ( Komossa &: Bad^ 19991 : Esauei et al. 20081) . For other TDEs, one 
did not detect the sources in both the rising and the declining phases, so that the peak luminosity cannot 
be well constrained. The maximum luminosity during the detected phase is only the lo wer limit of the pea k 


luminosity. However, since the lightcurve peak is rounded and spreads in several weeks (|Lodato et al 
the observed peak luminosity may not be too different from the true peak luminosity. 


i(^ . 


Due to the narrow bandpass of the X-ray telescopes that detect TDEs, the spectral shape for TDEs is 
not well constrained. Except for the BAT detected TDEs for which a treatment similar to GRBs can be 
applied, for the majority of TDEs, we only apply an empirical relation to estimate the bolometric luminosity 
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by multiplying X-ray luminosity by a factor of 1.4, i.e. (IMaksvm et al.ll201(3r) 

Lp^hol 


= 1.4 X Lp . 


(30) 


For Swift TDEs, a power-la w spectrum in the BAT band is reporte d for both events with a photon index 
of 1.8 and 1.6, respectively (IBurrows et al.ll201ll : ICenko et al.ll2012ll . Even though the spectrum of jetted 
TDEs is not known, we speculate that they have a non-thermal spectrum with Ep not far above the BAT 
band. As a result, we apply Ea. (l25l) to estimate Ep and apply Ea. (l2^ to estimate the fc-correction factor 
assuming the standard Band-function parameters for the spectrum. This led to a fc-correction factor 2.1, 
and 2.2 for Sw J1644-I-57 and Sw J2058-I-05, respectively. 


5. Results 
5.1. LGRBs 


Since LGRBs and SGRBs have different physical origins (massive star core collapse vs. compact star 
mergers), we derive their event rate density and luminosity function separately. For each type, we first 
derive the GLF by ignoring possible redshift evolution effect. Then we dedicate one subsection to discuss the 


in the luminosity function (Liang et al. 2007; Virgili et al. 

2 OO 9 I). which may have a somewhat different 

ohvsical origin fe.g. Liang et al. 2007; Bromberg et al. 2011; 

Nakar & Sari 20121. In our analvsis. we adoot 

two approaches. First, we derive the luminosity function 

jf HL-LGRBs and LL-LGRBs separately (e.g. 
b, we derive the luminosity function of LGRBs 
e. TPL) luminosity function. This approach is 

Liang et al. 2007; Virgili et al. 2009ll. In the second approac 

jointly by fitting the data together with a two-component (i 


justified in view that both LL- and HL-LGRBs are associated with Type Ic supernovae, and therefore may 
share a common physical origin. 


5.1.1. Global luminosity function 


Figure [3] shows the results of GRBs. The event rate density above a particular L as a function of the 
bolometric luminosity L is presented in the left panel, and the specific event rate density as a function of L, 
which describes the GLF, is presented in the right panel. In both panels, the HL-LGRBs (red), LL-LGRBs 
(blue), and SGRBs (black) are presented separately. The luminosity bin is taken as 0.3 (LGRBs) or 0.7 
(SGRBs) in the logarithmic space. The horizonta l errors denot e the luminosity bins, whereas the vertical 
errors are calculated from small-sample statistics ( Gehrels 1986[) . The best fit lines for all three sub-types 
of GRBs are also presented. The fitting results are summarized in Table [S] and Table [T] For all the fitting 
parameters, we also present the Icr range of errors based on 5000 sets of Monte Carlo simulations. 


The local event rate density of LL-LGRBs is 164lg® Gpc 
.,-1 


^ yr ^ 


with a minimum luminosi ty 5 x 


10'^®ergs“^, which is roughly consistent with previous results ( Liang et al. 2007llVirgili et al. 2009ll . The 
longer working period of Swift till now than the ones relevant for the previous two papers makes the event 
rate density slightly lower than before. However, it is still around two orders of magnitudes higher than that 


of HL-LGRBs 


which is 2.4+°;| Gpc-^yr"! 


above 3 x lO'^^ergs”^, or 0.8^°;} Gpc“^ yr“^ above 10®°ergs“^ 
(a typical luminos ity t hreshold adopted befo r e). T his local event rate of HL-LGRBs is slightly lower than 
Liang et al. ( 2007 1 and Wanderman fc Piran ( 2010l l. 
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From Eq. (USD, one can see that po,L is proportional to the luminosity function $(-F), with po defining 
the normalization. We use either a single power law (SPL) or a triple power law (TPL) to ht the po.L 
distribution for each sub-type of GRBs to characterize their GLFs. The htting results are shown in Table 
[SI HL-LGRBs can be ht with a TPL with = 2 . 2 lg' 2 , = l.Olg}, = 2 . 0 ;to '3 break 

luminosity L™ = 5 . 0 t ?;3 x 10®° erg s“^, L™ = T.llg;^ x 10®^ erg s“^. The luminosity distribution of LL- 
LGRBs can only be ht using a SPL with = 2 . 3 lg' 2 . The LL-LGRBs are not the straightforward extension 
of HL-LGRBs to low luminosities. But the slope of LL-LGRB is similar to that of hrst component of HL- 
LGRBs. The normalization of LL-LGRBs is a little bit lower as we apply a lower threshold for LL-LGRBs 
than HL-LGRBs. Therefore it may be possible that LL-LGRBs follow the extension of LF of HL-LGRBs 
if we take a TPL ht to the joint LL- and HL- samples. Such a ht is presented in FigjU with the best ht 
parameters being ai = 1.7^°'^, a 2 = I.OIq'ii cig = 2 . 0 ^°^ and the break luminosities Lb,i = l-O^oJ x 10 ®^ 
erg s“^, Lb ,2 = T.Sltg i x 10®^ erg s“^. This is similar to the TPL ht to HL-LGRBs alone except for a slightly 
shallower ai, which is compromised by the slight mis-match between LL-LGRBs and HL-LGRBs. 


The TPL nature of the joint G LF for LGRBs is interesting. Such a feature was only noticeable in the 
past when LL-LGRBs are included (jLiang et al.ll2007l : IVirgili et al.ll2009ll . With the current sample, we hnd 
that it is required even for HL-LGRBs alone. With the joint ht LL- and HL-LGRBs, we hnd that the steep 
component in the low-luminosity end now extends all the way to ^ 10 ®^ erg s“^, so that no clear separation 
between the LL- and HL-LGRB population is seen. 


We would like to stress that if we focus on the two high-L segments in the GLF, our results (q !2 1.0, 

as ~ 2.0, Lb, 2 7.8 x 1 0®^ erg s~^) are br oadly consistent with previous results: ^ 0.8, ^ 2.6, 

and ^ 2.5 x 10®^ bvlLiang et alJ (|2n07ll : and a^^ - 1.17, - 2.44, and L^^ - 3.1 x 10®^ ergs"" by 

Wanderman fc PiranI ( 2010[ ^ 


According to Eqs.([T|) and (fTTjl . the indices of po.>L and po.udL should be roughly the same, i.e. the 
index of po,>L should be greater than the index of po,L by one. This is generally satished for most of the 
transients studied in this paper (see all the indices marked in Figs.3-8). 


5.1.2. Luminosity function evolution 

Using method laid out in Section [221 we study the evolution effect of LGRB luminosity function. The 
results of pg l are shown in FigjSjleft panel. For each redshift bin, we fit the luminosity function with a SPL 
or BPL if the latter is needed. The fitting parameters are presented in Table [H One can see that indeed 
there is an apparent luminosity function evolution effect. However, there is no clear pattern to quantify the 
evolution. The right panel of Fig|5] shows the redshift-dependent break luminosity. In some redshift bins a 
break is clearly seen. However, in some other bins, the break either does not exist, or simply required by 
only one data point with low significance. For SPL fits, we place either a lower limit or an upper limit based 
on the highest or lowest luminosity bin. As shown in FigjSj right panel, there is no clear pattern to quantify 
the evolution effect. Further more, the PL indices ai and a 2 also show signihcant evolutions (variations) 
in different redshift bins. Since different redshift bins have different L^, and since the behavior below L^ 
is poorly constrained by the data, in Table | 8 | we choose different L^ for different redshift bins. For nearby 
universe, we also get a p§ at higher L^. We obtain a consistency of pg derived from data from different 


^Notice and in previous works correspond to 02 and Q3 in our notation. Also I Wanderman &: PiranI ll201Cll ) derived 
the luminosity function in the logarithmic space, so the a and 0 values in their notation are smaller by one from our values. 



























- 16 - 


redshift interval. 


To minimize the truncation effect by the flux limit of detectors, we also use a sub-sample with a higher 
threshold (with peak photon flux larger than 1.8 ph s“^ cm“^), see FigEl The sub-sample consists less 
GRBs in some redshift intervals (e.g. z S (0.5,1.5)) so that the luminosity function can be fit by a SPL. In 
the nearby universe, on the other hand, since LL-LGRBs are dropped due to this high-threshold criterion, 
the luminosity function demands a TPL fill^. The fitting parameters are also presented in Table |51 The 
evolution of Lf, is now better quantified by oc (1 -I- z)^ ’’ (FiglHlright panel), but a signficant variation of 
ai and a 2 values in different redshift bins remain. We therefore conclude that there is no straightforward 
way to quantify the the evolution effect of LGRB luminosity function. 


On the other hand, if we assume the evolution law assumed from the latest papers fe.g. lPetrosian et al 


2ni5c IYu et al.N2ni5c IPescalli et al.N2015[l . i.e. the PL indices remain unchanged, and only Lf, evolves with 


k ~ 2.3, we can map the luminosity function from the data by correcting the luminosity L at z to the “local” 
value, i.e. Lq = L/(l + z)^-^, to derive the “local” luminosity function. The result is shown in FiglTl This 
local luminosity function can be fit by a broken power law with ai — 1.5, a2 = 2.5 and Lb o = 51.6 erg s“^. 
This is consistent with the results presented in previous papers. 


5.2. SGRBs 


5.2.1. Global luminosity function 


luminosity 7x 10^®ergs it is 4. 7.lt2-2 Gpc- ■3 yj- 1 for the Gaussian, log-normal and power- 


law delay models, respectively. Taking a typical minimum luminosity 10®® erg s ^ 


The local event rate density for SGRBs vary slightly for different merger delay models. For a minimum 

yr" 

as adopted by previous 

authors, the numbers are 1. 3 ^ 0 . 3 ) l-2to 3 , and 3.3t|!|;g Gpc“®yr“^, re spectively. The local event r ate density 
for the log-normal model is slightly lower than the value reported by IWanderman fc PiranI (j2ni4l l , as we use 
a slightly lower flux threshold for Swift BAT-detected SGRBs. 

SGRBs come with a much smaller number of events than LGRBs. Assuming that all the SGRBs with 
redshift measurements are of a compact star merger origin, we derive their luminosity function in FigISl The 
three different merger delay models give slightly different results, but in general all three models are consistent 
with having a SPL luminosity function with an index of ^ 1.6 (for details, see Figlfll). This is different from 
HL-LGRBs which require a TPL luminosity function. It is also different from IWanderman fc PiranI (|2014ri , 
who claimed a BPL luminosity function. 

One caveat of our treatment is the as sumption that all SGRBs are of a compact star merger (Type I) 
origin. In the Swift era. IZhang et al.l (j2009l l first suspected that some high-L SGRBs at high redshifts may not 
be of the Type I origin, but may rather originate from massive star core collapse (Type II). They suggested 
to apply mu lti-wavelengt h crite ria (instead of using y-ray duration only) to judge the physical origin of a 
GRB. Later Virgili et al.l ( 2011 1 pointed out that the assumption that all SGRBs are of the compact star 
merger origin (Type I) is disfavored since one cannot simultaneously account for the z — L two-dimensional 
distribution and the log TV — logP distribut i on of SGRBs. They claimed that a good fraction of SGRBs 
may be of a Type II origin. Bromberg et al. ( 2013[) recent ly reached the similar conclu sion using a different 
argument based on the duration distribution of GRBs. IWanderman fc PiranI (|2014ll used the criteria of 


®For the full sample, the inclusion of LL-LGRBs compensates the low-L excess so that a BPL presents a reasonable fit. 
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Bromberg et al. ( 20131) and excluded about 1/3 of SGRBs that they suspect to be of a massive star core 


collase origin. This may explain the difference between the results of this paper (SPL luminosity function) 
and that paper (BPL). We did not exclude any SGRB from our sample for the following reason. We believe 
that any conclusion a bout the physical category of a SGRB based on the duration informati on only (e.g. 
Bromberg et al. 2013ll is not reliable. Rather one should consider multi-wavelength criteria (jzhang et al 


200911 . espe cially the host galaxy type an d the afterglow location within the host. Host galaxy studies 


of SGRBs ( Fong et al. 2010t Berger 2014 1 suggested that the hypothesis that all SGRBs belong to the 


compact star merger (Type I) category is not in conflict with the data. Indeed, some SGRBs excluded by 


Wanderman fc PiranI (1201411 (dehned by them as collapsars) actually haye large offsets from host galaxies, 
fully consistent with being due to a compact star m erger origin (e.g. GRB 070724 with offset 5.46 ± 0.14 


kpc, and GRB 070809 with offset 33.22 ± 2.71 kpc, iFong et al 


2010 ; 


Fong fc: Bergen 1201311 . Furthermore, 


a recent study of SGRB emission amplitude parameter (ILii et al.ll2014^ also suggested that most obseryed 
SGRBs are not disguised SGRBs due to the “tip-of-iceberg” effect. Eyen though we belieye that there exists 
a contamination Type II GRBs in the SGRB sample, without studying the multi-wayelength data of SGRBs 
in detail, we belieye that it is more appropriate not to exclude any SGRB in this study. 


5.2.2. Luminosity function evolution 

For short GRBs that haye much less data, it is more difficult to study the possible eyolution effect 
of luminosity function. We manage to diyide the short GRB sample into three redshift bins {z < 0.3, 
0.3 < z < 0.6, and z > 0.6), and apply the method in Section [2(^ to deriye pg using the data in different 
redshift bins. The results are shown in FiglU The luminosity functions can be all fit with a SPL. The slope 
in the first redshift bin (0 < z < 0.3), l.Oj^g'g, is somewhat steeper than the those of other two redshift bins 
(1-5 -o’2 0.3 < z < 0.6; and l.dlgj for z > 0.6, respectively). However, the slopes are consistent with 

each other within error. 


5.3. Shock breakouts 


The results for these two SBO events are presented in FigHOl In view of the possible connection between 
SBOs and LL-LGRBs in general, we present LL-LGRBs in the same plot for comparison. 

Our results suggest that the event rate density is 3 .lt.2.o ^ Gpc“^yr“^ for XRO 080109/SN 2008D- 
like SBO events (luminosity ^ 6.1 x 10^^ erg s“^), and is Gpc“^yr“^ for GRB 060218/SN 2006aj-like 
SBO events (luminosity ^ 1.5 x 10"'’^ erg s“^). The former implies that the local event rate density of SBO 
is at least 10^ times higher than that HL-LGRBs. A sensitive large field-of-view X-ray detector would lead 
to discovery of a large sample of these events. 


For both hgures we find that XRO 080109/SN 2008D (green data point at low luminosity) follows 
the extension of LL-LGRBs (blue line). We also perform a joint fit between SBOs and LL-LGRBs and 
get a SPL luminosity function with slope ~ 2.0, which is similar to that of the slope of LL- LGRBs only 


= 2.3). This lends support to the possible conn ection between LL-LGRBs and SBOs (e.g. IWang et al 


2007 : Nakar fc Sari 20121 : Barniol Duran et al. 2014 1. 
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5.4. Tidal disruption events 


Figure [TT] presents the results for TDEs. As mentioned above, two caveats are that the measured 
luminosities are only the lower limits of the “peak luminosities”, and that some uncertainties are associated 
with the fc-correction parameters. Bearing in mind these caveats, the following conclusions can be drawn. The 
event rate densities at different luminosity bins have a wide distribution, ranging from over 10® Gpc“® yr“^ 
at 10'^^ erg s“^ to ^ 10^ Gpc“® yr“^ at 10^® erg s“^. In the luminosity range of (10"^® — 10"^'*) ergs“®, 
the event ra te density is in the ran ge of (10® — 10"^) Gpc“® yr“®, which is co nsistent with both theoretical 


predictions (Wang fc Merritt 2004) and the estimates based on observations ( Esauei et al. 20081 : Luo et al 


2008; Maksvm et al. 2010ll . The event rate density of the two Swift TDEs is ^ 0.03i'^Q g 2 


Gpc“®yr“®. Similar 

to GRB 060218, these two events were de tected through image triggers. In particular, Sw J1644+57 triggered 
BAT multiple times (all image triggers) (IBurrows et al.ll201l[ l. We use the brightest peak to define the peak 
luminosity, and derive the event rate density based on the trigger information of that epoch. 


The GLE of TDEs (including both normal TDEs and jetted TDEs) can be roughly described by a 
single PL with = 2.0. With the current sample, the jetted TDEs detected by Swift seem to lie in the 

extension of normal TDE luminosity function to high-luminosities. The event rate density of Sw J1644-I-57 
shows a flattening feature at the highest luminosity. More data are needed to verify whether jetted TDEs 
form a distinct component in the global TDE luminosity function. 


5.5. Global distribution 

In Fig[T2j we plot all high-energy transients in one figure. The left panel shows the local event rate 
density above a certain luminosity {po,>l) as a function of L, and the right pabel presents the specific local 
event rate density (po.n), which essentially represents the GLF. All the data points are given a Icr error in 
vertical axis. The width of luminosity bin is shown as the horizontal error bar. 

Intriguingly, all the transients seem to line up to form a rough single power law distribution. A best fit 
to the GLFs of all the transients gives a slope of = 1.6. All the events lie within the 3 ct confidence 

bounds of the best fit. The region below this correlation line is likely due to an observational bias, and could 
be filled with new types of transients. With a much lower event rate density, these transients may not have 
been detected within the time span of the modern high-energy astronomy. The region above the correlation 
line, on the other hand, is not subject to selection effects and must be intrinsic. The existence of such an 
upper boundary of high-energy transients is intriguing, which may be rooted from more profound physical 
reasons. 


6. Conclusions and discussion 

In this paper, we systematically investigated the local event rate density, redshift evolution, and GLF 
of several known extra-galactic high-energy (y-ray and X-ray) transients, including high-luminosity and low- 
luminosity long GRBs that have a massive-star core-collapse origin, short GRBs that likely have a compact 
star merger origin, supernova shock breakouts, and tidal disruption events of stars by super-massive black 
holes. Our conclusions can be summarized as follows. 


• For all types of transients, the GLFs are typically well described by a single power law, although that 
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of HL-LGRBs demands a triple power law. 


• The local event rate density of each type of transient depends on the minimum luminosity. For GRBs, 
we get = 164;®® Gpc-®yr-i for = b x 10'^^ ergs"! for LL-LGRBs; = 2.4;°;® 

Gpc-®yr-i with = 3 x 10^® ergs"®) for HL-LGRBs, and = 4.2;^^, Tgfj;®, T.i;?;® 

Gpc“®yr“® for SGRBs, with = 7 x lO"*® ergs“® for the Gaussian, log-normal and power law 

merger delay models, respectively. Even though with two confirmed cases, the SBOs have event rate 
densities cover a wide range, from 3.1'^2'q ^ 10"^ Gpc“®yr“® for XRO 080109/SN 2008D-like events 
(luminosty ~ 6.1 xlO"'® erg s“^) to li;g® Gpc“®yr“^ for GRB 060218/SN 2006aj-like events (luminosity 
~ 1.5 X lO'^^ erg s“®). The event rate density of TDEs also covers a wide range, from l.O;®; x 10° 
Gpc“®yr“® with L™® = 10'^^ ergs“® for normal TDEs to 0.03 ;q'q 2 Gpc“®yr“® above 10"*® ergs“® for 
jetted TDEs detected by Swift. 


• For GRBs, we confirmed the previous work (|Liang et al.ll2007l : IVirgili et al.ll2009ll that LL-LGRBs still 
do not straightforwardly follow the extension from HL-LGRBs. However, a TPL fit to the entire LGRB 
population suggests that the steeper GLF slope in the low-energy end now extends to a much higher 
luminosity, so that LL- and HL-LGRBs are no longer clearly separated. Consider the GLF shape in 
the high-L end, the indices {a = 1.0, P = 2.0) and the break luminosity (Lb ,2 = 7.8 x 10°^ erg s~®) a re 


generally consistent with (even t hough not identical to) what was found in previous work (ILiang et al 


2007 : Wanderman fc PiranI 20in[) . For SGRBs, w e found a SPL luminos i ty fun ction with q-SGRB _ 


in contrast with the BPL distribution found bv IWanderman fc PiranI (120141) . The discrepancy may 
lie in different sample selection criteria: whereas IWanderman fc PiranI (1201411 excluded about 1/3 of 
SGRBs, we included all the SGRBs in our analysis. 


We c onfirm the conclusion of previous authors (e.g. Petrosian et al. 2015 : Yu et al. 20151: Pescalli et al 


20151 1 that the luminosity function of LGRBs likely evolve with redshift. However, we find that the 


evolution cannot be easily quantified with a simple analytical model. Nonetheless, if one assumes that 
shape of the LF does not change and only the break luminosity evolves with redshift, the consistent 
results as previous authors can be achieved. 


• SBOs and LL-LGRBs have a similar index = 2.0, supporting the idea that LL- 


LGRBs mav be related to shock breakouts (Camoana et al. 

20061 

Waxman et al. 

2007; 

Wansf et al. 

2007; 

Bromberff et al. 

2011 : 

Nakar & Sari 

20121 . 


• The global luminosity function of TDEs is consistent with a single power law with = 2.0. The 

jetted TDEs discovered by Swift seem to be consistent with the extension of normal TDEs to high- 
luminosty regime, even though a flattening feature is seen. More data are needed to judge whether 
jetted TDEs form a new component in the GLF. 


• Intriguingly, all the high-energy transients are consistent with having a global single power-law distri¬ 
bution of GLFs with a slope 1.6. Even though there could exist transients below the line which have 
not been discovered, the lack of events above the line is real. The existence of such an upper boundary 
is intriguing, and its physical origin is unknown. 


To perform this analysis, we adopted/derived the redshift distribution factor f{z) of various types of 
transients based on different models. For LGRBs and SB Os, we assume that the event rate density 
traces the SFH, and adopt the empirical model fEa. dTOl) ! of Yiiksel et al. ( 2008ll . For SGRBs, through 
Monte Carlo simulations, we derived empirical f{z) functions for three merger delay models: Eas. d^ . 
(ED, and (1^^ for the Gaussian, log-normal, and power-law, respectively. For TDEs, we assume that 
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the event rate is constant within each galaxy, and deri ved an empirical form ula of f{z) lEa. (l24l) l based 
on the black hole number density evolution following IShankar et al.l (12013^ . These empirical formulae 
can be directly used in the future. 


This paper focuses on high-energy transients only, serving as a reference for t he wide-held y-ra y as¬ 
tronomy and the upcoming wide-held X-ray astronomy (led by e.g. Einstein Probe. lYuan et al.ll2015ll . We 
notice that the phenomenology of all the transients studied in this paper extends to lower frequencies in the 
electromagnetic spectrum. For example, GRBs have multi-wavelength afterglows. SBOs can peak in the UV 
or even optical band if the progenitor star is large enough. UV and optical TDEs and radio counterparts 
of jetted TDEs have been discovered. The study of high-energy transients in the low-frequency domain is 
beyond the scope of this paper. 


We thank an anonymous referee for constructive comments and suggestions, Takanori Sakamoto for im¬ 
portant discussions on Swift trigger algorithms, Francisco J. Virgili for important communications on Monte 
Carlo simulations, Hou-Jun Lii for providing the data of Epeak of short GRBs, Bin-Bin Zhang for technical 
help on data processing, Junchen Wan for communications on statistics, George Rhee and Kayhan Gultekin 
for discussions on black hole number density redshift distribution, and Amy Lien for helpful communications. 
This work is partially supported by National Basic Research Program (973 Program) of China under Grant 
No. 2014CB845800, the NSEC (11273005) and SRFDP (20120001110064). Hui Sun is supported by China 
Scholarship Program to conduct research at UNLV. 
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Table 1: Instrumental parameters. 


Detectors 

Operation Time 

Field of View 

Sensitivity (erg cm ^s 4) 

(Instrument) 

(T) 

(0) 

(^’th) 

CGRO(BATSE) 

10 yrs 

TT sr. 

3.0 X 10-® 

HETE-II(WXM) 

7 yrs 

0.8 sr. 

8.0 X 10-4* 

INTEGRAL (IBIS) 

12 yrs 

0.26 sr. 

9.1 X 10-® 

Swift (XRT) 

10 yrs 

5 X 10“®sr. 

10-42 (1000s) 

ROSAT(PSPC) 

8 yrs 

10“3 sr. 

3.0 X 10-43 (500s) 

XMM-Newton(EPIC) 

15 yrs 

2 X 10-4 sr. 

2.0 X 10-44 (lO^s) 

Chandra(AGIS) 

14 yrs 

6 X 10-4 sr. 

4.0 X 10-45 (lO^s) 

Swift (BAT) 

10 yrs 

1.33 sr. 

3 X 10-® for HL-LGRBs and rate-triggered LL-GRBs 

Swift (BAT) 

10 yrs 

1.33 sr. 

lO-"^ for SGRBs 

Swift (BAT) 

10 yrs 

1.33 sr. 

2.8(3.1) X 10-® for LL-LGRB 060218/100316D 

Swift (BAT) 

10 yrs 

1.33 sr. 

10-® for Sw J1644-h57 

Swift (BAT) 

10 yrs 

1.33 sr. 

4.3 X 10-44 for Sw J2058-b05 


Table 2: Best-fit merger delay models of SGRBs with respect to star formation history. 


Delay Model 

Formula 

Best-fit parameters 

Reference 

Gaussian (G) 

mG{T)dT = exp j l\/^at,GdT 

tt,G = 2 Gyr, crt,G = 0.3 

(1) 

Log-normal (LN) 

TOLN(T)(ilnT = exp /(y^cr^_LN)dlnT 

tt.LN = 2.9 Gyr, CTt,LN = 0.2 

(2) 

Power law (PL) 

<7PL('r)dT = T ‘^^dr 

at=0.81 

(2) 


References. — I'll lVirgili et al.l ll201lh : f21 IWanderman fc Piraiil <l2014h . 
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Table 3: The LL-LGRB sample. 


Name 

Detector 

Energy band 

-Ppeak 

T h 

4^p,bol,48 

redshift 

TP c 

-^peak 

Reference‘s 

GRB 980425 

CGRO 

50-300 keV 

0.96 ±0.05 

0.058 

0.0085 

122 keV (GPL) 

(1) 

XRF 020903 

HETE II 

2-10 keV 

2.2 ±0.8 

7.42 

0.251 

2.6 keV (GPL) 

(2) 

GRB 031203 

INTEGRAL 

20-200 keV 

1.3 ±0.0 

9.85 

0.155 

121 keV (PL) 

(1) 

GRB 060218 

Swift 

15-150 keV 

0.25 ±0.11 

0.147 

0.033 

4.5 keV (GPL) 

(3) 

GRB 080517 

Swift 

15-150 keV 

0.6 ±0.2 

3.03 

0.09 

202 keV (PL) 

(4) 

GRB 100316D 

Swift 

15-150 keV 

0.1 ±0.0 

0.116 

0.0591 

19.6 keV (GPL) 

(5) 


“Peak photon flux in unit of phcm The values for GRB 980425 and GRB 031203 are taken from GCN 67 and 

GGN 2460 seperately. XRF 020903 is from the ISakamoto et all ll2004h . Swift samples are downloaded from swift table 
(http : !Iswift.gsjc.nasa.gov!archive!grhjtahlel). 

^Peak bolometric luminosity is calculated after k-correction. Tp^bol,48 unit of lO'^® ergs”^. 

is either dire ctly given from spectrum fit in literatures for Power-law with a cutoff (CPL)-fit models or calculated through 
IVirgili et alj |2^^ for Power-law (PL)-flt models. 

References. — ( l) fKaneko et al.l ll2007h : ('2i lSakamoto et al.l ll2004li : f3l lCampana et al.l ll2006lh Gl lStanwav et ^ flo^); 
fSi lPan et al.l l|2Qllh 


Table 4: The shock breakout sample. 


Name 

Detector 

Energy band 

T a 

-f^p,bol,46 

redshift 

Reference 

SN 2006aj/GRB 060218 
SN 2008D/XRO 080109 

Swift (BAT) 
Swift (XRT) 

15-150 keV 

2-10 keV 

14.7 

0.0061 

0.033 

0.007 

Gamnana et al. ('2006') 

Soderberff et al. ('2008') 


The bolometric luminosity of XRO 080109 is calculated from fc-correction based on a power-law spectral with photon index of 
2.3. Tp^bol,46 is in unit of 10^® ergs“^ 
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Table 5: The TDE sample. 


Name 

Detector 

Energy band 

T ^ 

-^p.bol 

redshift 

Reference 

NGC 5905 

ROSAT 

0.1-2.4 keV 

3.6 X 

10^3 

0.011 

Bade et al. ('1996'): Komossa & Bade (1999') 

RX J1242 

ROSAT 

0.1-2.4 keV 

1.2 X 

1044 

0.05 

Komossa & Greiner ('1999') 

RX J1624 

ROSAT 

0.1-2.4 keV 

2.2 X 

1044 

0.064 

Grime et al. ('1999') 

RX J1420 

ROSAT 

0.1-2.4 keV 

7.6 X 

10^3 

0.147 

Greiner et al. (20001 

NGG 3599 

XMM 

0.2-2.0 keV 

7.1 X 

lO^i 

0.0028 

EsQuei et al. (2007') Esauei et al. (2008') 

SDSS J1324 

XMM 

0.2-2.0 keV 

6.7 X 

10^3 

0.088 

Esauei et al. (2007') Esauei et al. (2008') 

TDXFJ 1347 

ROSAT 

0.3-2.4 keV 

8.8 X 

1042 

0.037 

GaoDelluti et al. (2009') 

SDSSJ 1311 

Ghandra 

0.3-3.0 keV 

7.0 X 

1042 

0.195 

Maksvm et al. (2010') 

2XMMiJ 1847 

XMM 

0.2-2.0 keV 

3.9 X 

1043 

0.035 

Lin et al. (2011') 

SDSSJ 1201 

XMM 

0.2-2.0 keV 

4.2 X 

1044 

0.146 

Saxton et al. (2012') 

WINGSJ 1348 

Ghandra 

0.2-2.0 keV 

2.8 X 

1042 

0.062 

Maksvm et al. (2013') 

Swift J1644-I-57 

Swift 

15-150 keV 

7.2 X 

1048 

0.354 

Bloom et al. (20111: Burrows et al. (20111 

Swift J2058-f05 

Swift 

15-150 keV 

7.6 X 

1047 

1.185 

Genko et al. (2012) 


“i'p.bol is in unit of ergs 


Table 6: The best fit luminoity functions of different types of extra-galactic high energy transients. For 
SGRBs, the results for three merger delay models (Gaussian (G), log-normal (LN) and power law (PL)) 
are given. The Ict errors of all the fitting parameters are presented based on 5000 sets of Monte Garlo 
simulations. 


Type 

Fit model 

ai 

Q!2 

03 

Lb (ergs 4) 

HL-LGRBs 

TPL 

o O+0.4 

^.z_o 2 

1-0+} 

2.0t°;3 

5.013;° 

X 10°°, 7.11+ X 1032 

LL-LGRBs 

SPL 

o O+0.2 

Z.0-0.2 

- 

- 


- 

Joint HL-/LL-LGRBS 

TPL 

1 '7 + 0.1 
-‘-•^-0.1 

l-O+i 

2.0+3 

l-O+.a 

X 1031, 7.81+ X 1032 

SGRBs(G) 

SPL 

-| f^+O-OS 
'‘-■'-0.08 

- 

- 


- 

SGRBs(LN) 

SPL 

-‘-■^-0.08 

- 

- 


- 

SGRBs(PL) 

SPL 

1 c+0.08 
'‘-■'-‘-0.08 

- 

- 


- 

Joint SBO/LL-LGRB 

SPL 

o n+0.09 
^•*^-0.09 

- 

- 


- 

TDEs 

SPL 

o n + 0.05 
^■^-0.05 

- 

- 


- 































-24- 


Table 7: The event rate density of various transients given an observed minimum luminosity threshold and 
a typical luminosity threshold. 


Type 

Lm (ergs"^) 

Po,>Lm (Gpc 3yr i) 

L' (ergs“i) 

Po,>L' (Gpc ^yr 

HL-LGRBs 

3 X 10^® 

2.4t0;i 

IqSO 


LL-LGRBs 

5 X 10^® 


1046 

440+??^ 

SGRBs(G) 

7 X 10^® 

4 2^^'^ 

^•^-1.0 

IqSO 

1 O+0.4 
-‘-•’^-0.3 

SGRBs(LN) 

7 X 10^9 

3 Q+1-2 

IqSO 

1 9+0.4 
-‘-•^-0.3 

SGRBs(PL) 

7 X 10^® 


IqSO 

o 0 + 1.0 

'^•’^-0.8 

SBOs 

1044 a 

X 

10^7 


TDEs 

10^2 

1.0t;^;3 X 10^ 

1044 

4.8t^;i X 102 

Swift TDEs 

1048 

n nQ+0.04 

U.UO_Q 02 

- 

- 


Since there are only two confirmed SBOs, two characteristic luminosities are given around the exact luminosities of the two 
events. 
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Table 8: The best fit Luminosity function parameters in different redshift bins for full-sample and sub-sample 
(Pp > 1.8 ph cm“^ s~^)- 


Full-Sample 


Redshift 

(X\ 

CX2 

Lb (erg s“i) 

P§.>L„ (Gpc-3 yr-i)- 

Lm (erg s“i) 

0 < z < 1 

1.6 

1.8 

4.0 X 10^^ 

1.6t°;3 (0.3t!];i) 

105O (1051) 

0.5 < z < 1.5 

1.4 

1.7 

^ 2.0 X 10^3^ oj. ^ 5 1 X 10^® 

i.2;°j (0.4tri) 

105O (1051) 

1 < z < 3 

1.1 

1.6 

1.2 X 10^2 

6 d 
+ 1 

cd 

1051 

2 < z < 4 

1.5 

> 4.0 

^ 5.5 X 10^3, or ^ 2.0 x lO’^^ 


1052 

z > 3 

~ 0 

1.7 

< 1.7 X 10^2 

0.2;0;i 

1053 


Sub-Sample 


Redshift 

(X\ 

a2 

Lb (erg s 1) 

Pg.>L„ (Gpc-3 yr-1) 

Lm (erg s 1) 

0 

< 


< 

1 

0.1'^ 

1.8 

8.5 

X 

1050 

i.7tr2 (o.Ttri) 

105O (1051) 

0.5 

< 


< 

1.5 

1.5 

- 

^ 1053, or 

^ 8.0 X 1050 

d d 
+ 1 

t- 

cd 

1051 

1 

< 

z 

< 

3 

1.1 

1.6 

3.2 

X 

1052 


1051 

2 

< 

z 

< 

4 

1.0 

1.5 

4.3 

X 

1052 


1052 


z 

> 

3 


~ 0 

1.9 

7.0 

X 

1052 

d d 
+ 1 

cd 

1053 


"^The local event rate density derived from the data in different redshift bins. The minimum luminosity varies at each bin due 
to the limited instrument sensitivity. We also give event rate density at same Lm for comparison. 

^In the redshift bin 0 < 2 < 1, LF could be fit by a TPL. Here and 0.2 denote the latter two components and Lf, is the second 
break luminosity. 
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Fig. 1.— The redshift distribution derived from Monte Carlo simulations for short GRBs considering three 
delay time models with respect to star formation history: Gaussian (black), log-normal (blue), and power- 
law (red). For each model, the result is derived from the average of 10000 simulations, each with simulated 
10000 systems. Dots are the simulated results, and the curve is the empirical multiple-power-law fits given 
in Eqs. dini), (HH) and (l22l) . 




Fig. 2.— Left panel: Redshift-dependent number density of supermassive black hol es with masses in th e 
range of 10® — lO^M© (magenta) and 10^ — 10® Mq (cyan) derived from the results of IShankar et al.l (|2ni3l l. 
Right panel: Normalized redshift distribution of TDEs through simulation (dots) and the empirical fit 
(Eq.dMl)). 
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Fig. 3.— Left panel: Event rate density (/Oo,>l) distribution for LL-LGRBs (blue), HL-LGRBs (red) and 
SGRBs (black). The luminosity bin has a width of 0.3 for HL-LGRBs, 0.5 for LL-LGRBs and 0.7 for 
SGRBs in the logrithmic space. For short GRBs, the Gaussian merg er delay time model is adopted. The 
vertical error bars represent the Itr Gaussian errors calculated from ( Gehrels IQSOll . The horizontal error 
bars show the width of the luminosity bin. Right panel: Luminosity functions of LL-LGRBs, HL-LGRBs, 
and SGRBs. LL-LGRBs and short GRBs can be fit with a single power law, with indices 2.3 and 1.7, 
respectively. HL-LGRBs are fit with a triple power law with oi = 2.2 = 1.0 and = 2.0. 



Fig. 4.— A joint fit of LL- and HL-LGRBs with a two-component luminosity function. 
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Fig. 5.— Left panel: Luminosity function of long GRBs in different redshift bins for the full sample. Blue 
dots denote the LL-LGRBs. Magenta, red, green, cyan and black dots represent HL-LGRBs from 0 < z < 1, 
0.5 < z < 1.5, l<z<3, 2<0<4 and z > 3, respectively. Best fit models are overplotted as solid curves 
with the corresponding color. Right panel: Break luminosity evolution inferred from the luminosity function 
fit from the left panel. The medium values of each redshift bin is taken. For dubious BPL, we give both 
upper limit (the minimum luminosity of data, lower triangle) and lower limit (derived from BPL fit, upper 
triangle). No clear evolution pattern is observed. 
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Fig. 6.— Same as FiglS]but for a sub-sample with peak photon flux larger than 1.8 ph s“^ cm“^. In the 
redshift bin 0.5-1.5 (red), the sample is fit by a single power law. The lower/upper limits of the break 
luminosity (red triangles in the right panel) are derived from the maximum and minimum luminosities of 
the data. These limits are not used in the L{, ht (the dashed line in right panel). 
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Fig. 7.— The derived local luminosity function of HL-LGRBs assuming a simple luminosity function evo¬ 
lution model with same luminosity function shape but an evolving break luminosity Lf, oc (1 -f The 

broken power law LF gives ai = 1.5, a 2 = 2.5 and Lf, = 51.6 erg s“^. 



Fig. 8.— Luminosity functions for SGRBs with three different merger delay models: Gaussian (black), 
log-normal (blue), and power-law (red). 
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Fig. 9.— Luminosity functions of short GRBs in three different redshift bins. 



log L {ergs ') log L (erg s ') 

Fig. 10.— Left panel; Event rate density (po.>l) distribution for SBOs (green) and LL-LGRBs (blue). 
Right panel: Joint luminosity function of SBOs (green) and LL-GRBs (blue). For both panels, single power 
law fits to LL-LGRBs alone (blue) are shown. One can see that the SBO event XRO 080109/SN 2008D 
roughly follows the extension of the blue line. For the right panel, we also show the SBO/LL-LGRB joint-fit 
luminosity function (red). One can see that the slopes of the blue and red lines are similar to each other. 
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Fig. 11.— Left panel: Event rate density (po,>l) distribution for TDEs. The luminosity bin has a width of 
1.0. Several results from previous works are also shown for comparison. Right panel: Luminosity function 
of TDEs with the best ht. 




Fig. 12.— Glocal distributions of all extra-galactic high-energy transients. Left panel: Event rate density 
above a minimum luminosity with respect to bolometric luminosity; Right panel: Joint luminosity function. 
Different types of events are marked in different colors. A single power law fit (red line) and 3cr boundary 
for the correlation are presented. 
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